{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Case study\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Yo-yo\n",
    "\n",
    "Suppose you are holding a yo-yo with a length of string wound around its axle, and you drop it while holding the end of the string stationary.  As gravity accelerates the yo-yo downward, tension in the string exerts a force upward.  Since this force acts on a point offset from the center of mass, it exerts a torque that causes the yo-yo to spin.\n",
    "\n",
    "![](diagrams/yoyo.png)\n",
    "\n",
    "This figure shows the forces on the yo-yo and the resulting torque.  The outer shaded area shows the body of the yo-yo.  The inner shaded area shows the rolled up string, the radius of which changes as the yo-yo unrolls.\n",
    "\n",
    "In this model, we can't figure out the linear and angular acceleration independently; we have to solve a system of equations:\n",
    "\n",
    "$\\sum F = m a $\n",
    "\n",
    "$\\sum \\tau = I \\alpha$\n",
    "\n",
    "where the summations indicate that we are adding up forces and torques.\n",
    "\n",
    "As in the previous examples, linear and angular velocity are related because of the way the string unrolls:\n",
    "\n",
    "$\\frac{dy}{dt} = -r \\frac{d \\theta}{dt} $\n",
    "\n",
    "In this example, the linear and angular accelerations have opposite sign.  As the yo-yo rotates counter-clockwise, $\\theta$ increases and $y$, which is the length of the rolled part of the string, decreases.\n",
    "\n",
    "Taking the derivative of both sides yields a similar relationship between linear and angular acceleration:\n",
    "\n",
    "$\\frac{d^2 y}{dt^2} = -r \\frac{d^2 \\theta}{dt^2} $\n",
    "\n",
    "Which we can write more concisely:\n",
    "\n",
    "$ a = -r \\alpha $\n",
    "\n",
    "This relationship is not a general law of nature; it is specific to scenarios like this where there is rolling without stretching or slipping.\n",
    "\n",
    "Because of the way we've set up the problem, $y$ actually has two meanings: it represents the length of the rolled string and the height of the yo-yo, which decreases as the yo-yo falls.  Similarly, $a$ represents acceleration in the length of the rolled string and the height of the yo-yo.\n",
    "\n",
    "We can compute the acceleration of the yo-yo by adding up the linear forces:\n",
    "\n",
    "$\\sum F = T - mg = ma $\n",
    "\n",
    "Where $T$ is positive because the tension force points up, and $mg$ is negative because gravity points down.\n",
    "\n",
    "Because gravity acts on the center of mass, it creates no torque, so the only torque is due to tension:\n",
    "\n",
    "$\\sum \\tau = T r = I \\alpha $\n",
    "\n",
    "Positive (upward) tension yields positive (counter-clockwise) angular acceleration.\n",
    "\n",
    "Now we have three equations in three unknowns, $T$, $a$, and $\\alpha$, with $I$, $m$, $g$, and $r$ as known parameters.  It is simple enough to solve these equations by hand, but we can also get SymPy to do it for us.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import init_printing, symbols, Eq, solve\n",
    "\n",
    "init_printing()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "T, a, alpha, I, m, g, r = symbols('T a alpha I m g r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "eq1 = Eq(a, -r * alpha)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "eq2 = Eq(T - m * g, m * a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "eq3 = Eq(T * r, I * alpha)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "soln = solve([eq1, eq2, eq3], [T, a, alpha])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "soln[T]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "soln[a]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "soln[alpha]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "The results are\n",
    "\n",
    "$T      = m g I / I^*   $\n",
    "\n",
    "$a      = -m g r^2 / I^* $\n",
    "\n",
    "$\\alpha = m g r / I^*    $\n",
    "\n",
    "where $I^*$ is the augmented moment of inertia, $I + m r^2$.\n",
    "\n",
    "You can also see [the derivation of these equations in this video](https://www.youtube.com/watch?v=chC7xVDKl4Q).\n",
    "\n",
    "To simulate the system, we don't really need $T$; we can plug $a$ and $\\alpha$ directly into the slope function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "radian = UNITS.radian\n",
    "m = UNITS.meter\n",
    "s = UNITS.second\n",
    "kg = UNITS.kilogram\n",
    "N = UNITS.newton"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:**  Simulate the descent of a yo-yo.  How long does it take to reach the end of the string?\n",
    "\n",
    "I provide a `Params` object with the system parameters:\n",
    "\n",
    "* `Rmin` is the radius of the axle.  `Rmax` is the radius of the axle plus rolled string.\n",
    "\n",
    "* `Rout` is the radius of the yo-yo body.  `mass` is the total mass of the yo-yo, ignoring the string.  \n",
    "\n",
    "* `L` is the length of the string.\n",
    "\n",
    "* `g` is the acceleration of gravity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = Params(Rmin = 8e-3 * m,\n",
    "                Rmax = 16e-3 * m,\n",
    "                Rout = 35e-3 * m,\n",
    "                mass = 50e-3 * kg,\n",
    "                L = 1 * m,\n",
    "                g = 9.8 * m / s**2,\n",
    "                t_end = 1 * s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a `make_system` function that computes `I` and `k` based on the system parameters.\n",
    "\n",
    "I estimated `I` by modeling the yo-yo as a solid cylinder with uniform density ([see here](https://en.wikipedia.org/wiki/List_of_moments_of_inertia)).\n",
    "\n",
    "In reality, the distribution of weight in a yo-yo is often designed to achieve desired effects.  But we'll keep it simple."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_system(params):\n",
    "    \"\"\"Make a system object.\n",
    "    \n",
    "    params: Params with Rmin, Rmax, Rout, \n",
    "                              mass, L, g, t_end\n",
    "    \n",
    "    returns: System with init, k, Rmin, Rmax, mass,\n",
    "                         I, g, ts\n",
    "    \"\"\"\n",
    "    L, mass = params.L, params.mass\n",
    "    Rout, Rmax, Rmin = params.Rout, params.Rmax, params.Rmin \n",
    "    \n",
    "    init = State(theta = 0 * radian,\n",
    "                 omega = 0 * radian/s,\n",
    "                 y = L,\n",
    "                 v = 0 * m / s)\n",
    "    \n",
    "    I = mass * Rout**2 / 2\n",
    "    k = (Rmax**2 - Rmin**2) / 2 / L / radian    \n",
    "    \n",
    "    return System(params, init=init, I=I, k=k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Testing `make_system`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "system = make_system(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "system.init"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write a slope function for this system, using these results from the book:\n",
    "\n",
    "$ r = \\sqrt{2 k y + R_{min}^2} $ \n",
    "\n",
    "$ T      = m g I / I^*  $\n",
    "\n",
    "$ a      = -m g r^2 / I^* $\n",
    "\n",
    "$ \\alpha  = m g r / I^*  $\n",
    "\n",
    "where $I^*$ is the augmented moment of inertia, $I + m r^2$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test your slope function with the initial paramss."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write an event function that will stop the simulation when `y` is 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test your event function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then run the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the final state.  If things have gone according to plan, the final value of `y` should be close to 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`theta` should increase and accelerate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_theta(results):\n",
    "    plot(results.theta, color='C0', label='theta')\n",
    "    decorate(xlabel='Time (s)',\n",
    "             ylabel='Angle (rad)')\n",
    "plot_theta(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`y` should decrease and accelerate down."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_y(results):\n",
    "    plot(results.y, color='C1', label='y')\n",
    "\n",
    "    decorate(xlabel='Time (s)',\n",
    "             ylabel='Length (m)')\n",
    "    \n",
    "plot_y(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot velocity as a function of time; is the yo-yo accelerating?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `gradient` to estimate the derivative of `v`.  How does the acceleration of the yo-yo compare to `g`?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
